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Abstract 

This note works out an advection-difTusion approximation to the den- 
sity of a population of E. coli bacteria undergoing chemotaxis in a one- 
dimensional space. Simulations show the high quality of predictions under 
a shallow-gradient regime. 
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1 Introduction 



I The chemotactic behavior of E. coli has been studied widely at both the mi- 

^ croscopic and macroscopic levels. The movement of bacteria involves a directed 

Jlf^ movement (run) and a random turning (tumble). Each individual carries an 

internal state which, in the presence of a time and space dependent external sig- 

psj nal, may be modeled as evolving according to a system of ordinary differential 

equations. In the presence of a signal (typically, a nutrient) in the environment, 

the individual changes its directiorTH at random, with a tumbling rate which de- 

_ pends on the internal state, biasing moves toward more favorable environments or 

^ away from noxious substances. This random reorientation introduces a stochastic 

• ^ character to the evolution equations, and the population behavior of such hybrid 

^ systems is modeled by jump-Markov state-dependent systems. 

H 

^ The transport equation that describes a jump-mar kov system is very difficult 

to study mathematically, and cannot be validated by typical experimental tech- 
niques such as optical density measurements of bacteria in microfluidics cham- 
bers. Thus it is of great interest to derive a simpler macroscopic equation for the 
density of bacteria from microscopic equations. In addition, the current interest 
in scale invariant transient behavior ("fold-change detection," see for example 
[3 m E]) requires such macroscopic descriptions when starting from the jump- 
Markov model, as discussed in [5j. Our goal in this note is to work out, in a 
shallow-gradient regime, the tools developed by Erban, Othmer, and Grunbaum 



^Potentially, its speed may change too, though this seems not to be a substantial factor in 
E. coli chemotaxis. 
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[U [2] for a mechanistically realistic model of E. coli chemotaxis. In the case 
of exponential gradients, the result is a constant-coefficient advection-diffusion 
equation. We provide the calculations as well as an agent-based simulation that 
verifies the theoretical predictions, for one-dimensional motions. Future work will 
expand these considerations to two and three dimensions. 



1.1 Preliminaries 



Let p{x, y, z/, t) be a density function describing a population of agents (for ex- 
ample, bacteria), modeled in a 2N + m dimensional phase space, where at time 

; we soon specialize to = 1) denotes the 
position of a cell, y = {yi, . . . , ym) ^Yc M™" denotes the internal dynamics of 
a cell, and v €V <Z denotes its velocity. Also, S{x, t) = {Si, . . . , Sm) £ M*^ 
denotes the concentration of extracellular signals in the environment. 

We assume that the following system of ordinary differential equations describes 
the evolution of the intracellular state, in the presence of the extracellular signal 
S: 

where /: x — t- is a continuously differentiable function with respect 
to each component, i.e., / G C^{W^ x M*^). 

The evolution of p with turning rate A = A(?/) is governed by the following 
transport (or "Fokker-Planck" or "forward Kolmogorov") equation: 

dp f 

-7^ + Vx ■ i^p + Vy ■ fp = -X{y)p + Ky)T{y, V, v')p{x, y, v\ t) du' (2) 

where the nonnegative kernel T(y, v, v') is the probability that the bacteria changes 
the velocity from v' to if a change of direction occurs. Also 



/ T{y,v,v') du = 1. 
Jv 



The main goal of this note is to derive a macroscopic model for chemotaxis using 
the microscopic model ([2]), i.e., we want to find an equation to describe the 
evolution of the marginal density: 



n(x,t) =11 p{x,y,i/,t) dydu. (3) 
Jv Jy 

As remarked in the introduction, this is of interest both because of experimental 
and theoretical reasons, in particular in the context of scale-invariant sensing [S]. 



2 Chemotaxis equation in one dimensional movement 

In this section, for simplicity, we study the movement of agents in one dimension 
and assuming a one-dimensional state, i.e., m = N = M = 1. We also take the 
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speed, u, as constant. Let p^{x,y,t) denote the density of the particles that at 
time t, are located at position x, with the internal state y, and with the constant 
speed v, and moving to the right (+) or left (— ) respectively. We assume that 
the following decay condition: 

p^{x,y,t) < A(x,^)e-"(-■'*)^ (4) 

for some functions A, a: M x [0, oo) — )■ M>o. The internal state evolves according 
to the following ODE system: 

^ = f^iy,,.,S,S'), (5) 

where /^:]Rx]Rx]Rx]R— s-M are continuously differentiable functions in each 
argument that describe the evolution of internal state of bacteria which move to 
the right (+) and left (— ) respectively. 

Note that we are allowing / to depend on the direction of movement as well as u 
and 5', the derivative of S with respect to space. In our example, = /~ only 
depends on y and S, but we can consider the more general dependence in these 
preliminary derivations. 

We consider the following equation for the tumbling rate: 

X{y,S,S') =giy,S,S'), (6) 
for some continuous function g. 



Then, according to Equation Q, p^{x,y,t) satisfy the following coupled first- 
order partial differential equations: 

^ + ^^ + = g{y,S,S'){-p^+p-) (7) 

+ §y[riy,'^^S,S')p~] = g{y,S,S')ip^-p-). (8) 

We now state the following lemma from [1] regarding the existence and uniqueness 
of solutions of 0-([8]): 

Lemma 1 (Existence and uniqueness of solutions). Suppose that G C"'^(M x 
R X R X R), and let S:M. x [0, oo) — )• R 6e continuous. In addition, assume 
that g in M) is always nonnegative, and : R — )• [0, oo) are given nonnegative 
compactly supported functions. Then there exists a domain Q C R x [0, oo) 
containing the entire line t = such that the system of equations {^-(^ with 
initial conditions p^ has a unique solution in Q. Moreover, the functions p^ 
are nonnegative wherever they are defined. 

The objective is to derive an equation for the macroscopic density function 

n{x,t) = / p'^{x,y,t) +p~{x,y,t) dy, (9) 
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using the microscopic model Q-Q, by the following technique from [T]. To this 
end we define the following additional moments: 



ni{x,t) = / {p^{x,y,t) +p {x,y,t)) dy, for i = l,2, ... 
j{x,t) = / u {p^{x,y,t) - p~{x,y,t)) dy, (10) 

..,../.v,....-..,.„.. ......... 



(11) 



Note that by condition Q all the moments are well defined. 
Next, we assume 

/+ = /o + ^^/l, 
/" = /o - ^fi, 
where the Taylor expansions of /o and fi are given as follows: 

fo = Ao + Aiy + A2y^ + ..., (12) 
/i = Bo + Biy + B2y^ + ..., (13) 

for some Ai's and i?j's that are functions of S, S' , and u'^. 
Also we consider the following Taylor expansion for g{y, S, S'): 

g{y,S,S') = ao + aiy + a2y'^ + (14) 

where a,'s are functions of S, S' . 

In addition, we assume = 0. Then by multiplying ([t]) and ^ by 1, u, and/or 
y, adding or subtracting, and integrating with respect to y on M, and applying 
the fundamental theorem of calculus and integration by parts, we obtain the 
following equations for macroscopic density and flux and their first moments: 

k>2 

^ + ^ = 5oi + Aim + Biji + + Bkjk, (17) 

k>2 k>2 

% + ^^^ = ^^Bon + u^Bim + (Ai - 2ao)ii (18) 



+ ^ BkUk + X](^fc - 2ak-i)jk 



k>2 k>2 

Note that by condition Q, for any i = 0,1, . . . we have that: 

lim y\p^ ±P~) = 0. 

y— >itoo 
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2.1 Parabolic scaling 



In this section, we introduce a parabolic scaling to derive a chemotaxis equation 



from the moment equations (15)-(18). Let L, T, vq, and A'^o be scale factors 
for the length, time, velocity, and particle density respectively, and define the 
following dimensionless parameters: 



n 



V 

n 



-- AiT, 



rii 



Bi 



3i 



rij 
No' 
BjL, for 



Ji 



i = 0,L 



for 



1,2, 



The parabolic scales of space and time are given by: 



eL \ X 



t 



,t_ 



(19) 

(20) 
(21) 

(22) 



for any arbitrary e. 

Now assume that under some conditions, for any i > 2, the jVs and n^'s are much 
smaller than ji and ni and can be neglected. (For example see Lemma [s] below.) 
Therefore, the dimensionless form of moment equations (15)-(18), become: 



^dn dj 

e — - + Cttt 
dt dx 



, dh 
dt dx 

^dfii dji 

e — - + e-— 
dt dx 

^dji ^2 "9^1 
e — - + -WT- 
dt dx 



0, 

-2aoi - 2aiji, 

eBoj + Aihi + eBiji, 

eu'^Bon + eu^Bifii + (^i — 2ao)ii 



(23) 
(24) 
(25) 
(26) 



Next, we write Equations (23)-(26) in a matrix form, as follows: 



,dw 



dt 



d_ 
' dx^ 



e" —r^ + e—Pw = eQw + Rw, 



(27) 



where w 



h, j ,ni, ji] and the matrices P, Q, and R defined as follows: 








1 





0\ 




1)2 






















1 


v 








i>2 


/ 



Q 

















\ 























Bo 





Bi 




v 


O'Bo 





u^Bi 





/ 



R 






Vo 












2ao 





-2a 





il 
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Assuming the regular perturbation expansion for w, 

where 



w 



n\f,n\,jlj , 



and comparing the terms of equal order in e in (27), we get: 

> = (n°, 0,0,0)"^ 



.0 



Rw 











(28) 



Rw'^ + Qw 
( 



7,0 



d_ 
dx 



Pw 



7,0 



-2aof - 2aii\ 
Ain\ 



Z>2 

ax 






(29) 



From Equation ( 29 ) , we can derive the following equation for j\ : 



jI 



and therefore, using the same Equation, we obtain the following equation for 
and j^: 



- 2aoi^ + 2ai-^ —h^ 

Ai - 2ao 

Now we compare the terms with order e^: 



^2 9 ^0 
ox 



Rw^ 



-Qw' + ^Pw' + ^.w'. 
ox dt 



(30) 



(31) 



Note that (1,0,0,0)"^ is in the kernel of R and the right hand side of (31) is in 
the image of R. Therefore their inner product is zero: 



9-1 (9 
ox dt 



0. 



(32) 



Equation ( 30 ) together with Equation ( 32 ) give the following equation for nP in 



the dimensionless variables: 



di 2ao dx? ao(ii - 2ao) 



(33) 



Since n{x,t) = n^{x,t) + 0{e), if we neglect the 0{e) term, Equation (33) leads 
to the following chemotaxis equation in dimensionless variables: 



dn 



92 



n 



aiBQu"^ dh 



di 2ao dx^ ao{Ai - 2ao) dx ' 



(34) 



Changing back to the original (dimensional) variables, we obtain the following 
PDE: 

dt ~ 2aodx^ ao{Ai - 2ao) dx' ^ ^ 
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3 Example 



In this example, we assume the internal state evolves according to the following 
ODE system: 

^ = PyiQ-a), (36) 

where a = Fr-rr, and p, o, K, and N are positive constants. 

l + i^(fr 

This system provides a simple model of chemotactic behavior in E. coli bacteria, 
as discussed below in the section on simulations. 

By ignoring the tumbling time, we consider the following equation for the tum- 
bling rate: 

A(y) = ra^, (37) 
where r and H are positive constants. 

The objective is to derive a parabolic equation for the macroscopic density func- 
tion. 

It is convenient to define a new internal state variable as follows: 

z = p{a-q). (38) 

A simple calculation shows that 

dz N . . N S' , , , , 

-r- = — z{z+pq){z+pq-p)±v w{z + pq){z + pq - p) (39) 

at p p b 

\{z) = ^{z+pq)^ (40) 
Lemma 2 (shallow condition). Let c = m.m.{pq,p — pq}. If 



< - and \z{Qi)\ < c. 



then \z{t)\ < c for all t >0. 



Proof. At 2;(0) = c, since 



5' c S' 

— < —, z ± V— > 0. By the definition of c, 
S S 

dz 

z + pq — p < 0, while z + pq > 0. Therefore at z = c, — < 0. On the other hand, 

dt 

S' 

at z{0) = — c, z ± u— < 0, z + pq — p < 0, and z + pq > 0. Therefore at z = — c, 

dz 

— > 0. Hence for any t > 0, \z(t)\ < c. □ 
dt 

We'll show that under the shallow condition, Lemma[2| and the following parabolic 
dimensionless parameters, the higher macroscopic moments j2, J3, • • •, and n2, ns, . . 
can be ignored. 
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Lemma 3. As before, let L, T, vq, and Nq be scale factors for the length, time, 
velocity, and particle density respectively, and define the following dimensionless 
quantities: 



e S' 



N = TN, p = p q = q and f = Tr. 



All other parameters remain the same as in Equations (19)- (20), and Equation 
(22). Then under the condition of Lemma^ for any i>l, 



— < ue and — < V^e , 
n n 



for some constants Ci = 0(1), andDi = 0(1). 



Proof. Note that 



< K implies 



< —K. Hence K 
e 



-K, where 



K = 0(1). Now by Lemma 2 1 we have \z\ < v — K = tvK. 



3i 



UqNo 
V 

V 



z^{v — V ) dz 



< — / z'{p'^+p ) dz {=vhi) 

< ^{eKvy I {p+ + p~) dz 



— Tt' C ' 

where Ci = K*i}'+'^ = 0(1). Note that — < A, where = — 

n V 



□ 



Using the notations of Equations (|12|)-(|13|) 



0, Ax = Npq (q-l), Bo = N-pq {q - 1) 



(41) 



and by notation of Equation ( 14 ) , the first two coefficients of the Taylor expansion 
of g is as follows: 



H rHq 

ao = rq , oi = 

pq 



H 



(42) 



Therefore, using the Equation (35), the chemotaxis equation for this particular 
example is: 



~di ~ 2^\q) dx^ ~ Npq (q-l)- 2rqH^dx' 



(43) 
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4 Simulations 



In this section, we provide agent-based simulations of the full jump-Markov sys- 
tem, and comparisons with the parabolic model, for systems of the special form 



in (36), which we repeat here for convenience: 

py[q ' 



The jump (or "tumbling" for bacteria) rate has the form X{y) = ra^ in (37). 
The parameters p, q, K, N, r, and H are all positive. 

For hgand concentrations Ki<^S<^Ka, where Kj ^ 18.2^M and Ka ~ 3000/iM 
are the dissociation constants for inactive and active Tar receptors respectively, 
the above equations provide a simple but phenomenologically accurate modej^ 
of the chemotactic response of E. coli bacteria to Me Asp; see for example [8], 
[3]. Furthermore, this is the range in which [7\ predicted, and [3j experimentally 
verified, scale-invariant behavior for E. coli responses to MeAsp. (See also [B] for 
further theory of scale invariance.) To stay in this range, we use ligand concentra- 
tions very close to, and mostly larger than, 5 = 100//M. Since our objective is to 
understand the quality of the parabolic (reaction-diffusion) equation, we depart 
slightly from models cited above, in postulating an instantaneous re-orientation 
after tumbling. A model with tumbling would require additional analysis. Also, 
since motion is one- dimensional, we ignore rotational random drifts from linear 
movement. 

The parameters in previous studies, see for example are as follows 

1/ = 0.0165, = 6, q = 0.5, r = 1280, H = 10, K = 0.000740, p = 0.017 

(in appropriate units corresponding to /iM concentrations, times in seconds, and 
lengths in millimeters). We start with these, but we will vary p in order to 
understand how the speed of adaptation (i.e, the time-scale at which the state 
variable y evolves) affects the quality of our theoretical predictions. 

In our numerical experiments, we take a one-dimensional channel of length 10 (in 
units of millimeters), and start all agents (cells) in the middle position, a; = 5, 
randomizing the initial direction of movement as right or left with probability 
1/2. The initial level y(0) of every cell is picked such that the activity a equals 
the adapted value q. (This is in accordance with the pre-adaptation setup in the 
microfluidics experiments in ^.) 

The length is picked large enough so that, in the time intervals considered (up 
to [0,200]), no boundaries are reached, so that, for all practical purposes, we 
are working on an infinite domain. We always take 100,000 cells, and display 
histograms based on 100 equal-sized bins. (These numbers represent a heuristic 
compromise between computational effort and smoothness of empirical densities.) 



^For convenience of analysis, we are using y = e"™ as a state variable, instead of the methy- 
lation level "m" as done in other papers. 

^In terms of the parameters used in [3], p — a{kR + kB), where a = 1.7 and kn = ks ~ 0.005. 
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We only consider exponential gradients S{x) = Ke^^, in which case the advection 
term does not depend on x, because S'{x)/S{x) = p is constant. We call p 
the "slope" (more precisely, this is the slope of log 5), and the shallow-gradient 
condition amounts to requiring p <^ 1. We typically pick p in the range 0.1 to 



1. Under the assumption that S{x) is exponential, our solution (43) becomes a 
constant coefficient advection-diffusion equation: 



dn 
dt 



D 



dx'^ 



V 



dn 
dx 



The general solution of this equation, when starting from a Dirac delta function 
at position xq{= 5), has the form 

n(x, t) = N{x + XQ- Vt, 2Dt) 

where M{p, a) is a Gaussian density with mean p and variance a. In other words, 
the solution is a translate of the fundamental solution of the heat equation. For 
purposes of comparison, the densities displayed at any given time t are plotted 
together with this theoretical prediction. As remarked earlier, the constant k is 
picked so that kS{xo) = 100, that is, k = 100e~^^. 



4.1 Slope p = 0.1, various values of p 

In Figures [T] to [5j we display simulated and theoretical distributions at times 
t = 10, 100, 200 as well as a plot of the means of the distribution on the in- 
terval [0,200]. Agreement to theory is very good, with larger p (faster internal 
adaptation dynamics) leading to closer fits. 






Figure 1: p = 0.1,p = 0.017 
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Figure 2: p = 0.1, p = 0.05 




Figure 3: p = 0.1,p = 0.1 
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Figure 4: p = 0.1,p = 0.5 




Figure 5: p = 0.1, p = 1 
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4.2 Slope p = 1, various values of p 

In Figures [6] to [9j we give plots with simulated and theoretical distributions at 
times t = 10, 100, 200 as well as a plot of the means of the distribution on the 
interval [0, 200]. Agreement to theory is now poor when p = 0.017 and p = 0.1, 
but is considerably better with larger p (faster internal adaptation dynamics). 




Figure Q: p = l,p = 0.017 

Observe that for low values of p, a large number of agents (cells) appear to be 
at the same "forward" position, moving to the right. An explanation of this 
phenomenon is that the probability of tumbling in each agent is very low. (The 
second peak is explained by the fact that 1 /2 of the cells where randomized to 
starting in a leftward motion. Thus, it takes a certain time for these cells to 
tumble and start moving right, toward higher nutrient concentrations.) 
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Figure 7: p = l,p = 0.1 




Figure 8: p = l,p = 0.5 
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4.3 Slopes p = 0.1 to 0.5, j9 = 0.017 



To further understand the behavior when p = 0.017, which matches theory well 
when p = 0.1, but badly when p = 1, we show in Figures 10 to 14 similar 
graphs for values p = 0.1,0.2,0.3,0.4,0.5. We display distributions at times 
t = 20, 30, 40. Asymmetry becomes more obvious for larger slope, and a small 
"right-moving front" can be seen arising at p = 0.5. 






Figure 10: p = 0.1, p = 0.017 



4.4 Slope p = 1, p = 0.017, 0.1, 1: plots of jump rates 



As remarked above, one may expect that the "right-moving front" observed for 
large slopes p and small adaptation time p is due to the jump (tumbling) rate 
A(t) (that is, X{t)dt the probability of a jump in an interval [t,t + 6t]) being low 
when p is small. This is indeed seen in the simulations. Figures [Tsj [TBI 17 
show the mean value of the jump rate A(t), averaged over all 100,000 cells in the 
simulation. Observe that the value approaches zero when p is small. However, 
for larger p, for example p = 1, these probabilities rapidly approach a more or 
else constant (and larger) value. These simulations are performed on the interval 
[0,100]. 

The mean and standard deviations of A(IOO) are, respectively, as follows: 



0.017 : 0.1286,0.7404; 0.1 : 0.4930,0.7413; 1 : 1.2232,0.2499. 
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Figure 11: p = 0.2,p = 0.017 




Figure 12: p = 0.3, p = 0.017 



17 



Figure 13: p = 0.4, p 



= 0.017 




Figure 14: p = 0.5, p = 0.017 
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viesri probability of turobling 




Figure 15: p = 1, p = 0.017, population means of A(t) 



viesri probability of turobling 




Figure 16: p = 1, p = 0.1, population means of A(t) 



mean probability of turobling 




Figure 17: p = 1, p = 1, population means of A(t) 
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4.5 Ideal steady-state values of A 



When S' / S = p, the z = p{a — q) variable in our theoretical derivation evolves 
according to the following cubic differential equation: 

dz N 

— = —{z+pq){z+pq-p){z±up) (44) 
at p 

where the + sign is picked when the agents are moving rightward and the — 
sign is picked otherwise. Physically, one is interested in solutions with positive 
so that the activity a{t) is always in the interval (0,1), which means, in 
terms of the z variable, that we must study Equation ( |44[ ) on the interval J = 
{—pq, —pq+p)- We will always assume that < q < 1, since otherwise the system 
has no equilibria for constant inputs (and in particular, does not perfectly adapt 
to step signals). Observe that J is forward- invariant, since i = when z = —pq 
and when z = p — pq. There is a third root of the cubic at z = ^vp, and this 
root belongs to the interval J if and only \i pq — p < zizi/p < pq. 

Suppose now that p > (source of nutrient is to the right), as in our simulations, 
and consider an agent that is moving rightward ("+" sign). We consider three 
cases: 

• PiQ — 1) < z^p < pq. in this case, as long as there are no jumps (direction 
reversals), z{t) z = —vp as t — )• oo. 

• pq < up: in this case, as long as there are no jumps, z{t) — t- —pq as t ^ oo. 

• up < p{q — 1): this case cannot happen, because q < 1 and up > 0. 

Thus, there is a bifurcation when the parameters satisfy: 

pq = up . 

In terms of the exponential rate for jumps X{t) 
first case that 

X{t) ^ A = -^{pq 
as t — )• cxD, and in the second case that 

\{t) -f . 

Now, suppose that an individual agent has spent enough time moving to the right 
that z has achieved a value close to its steady state. If the parameters are such 
that pq < up or if pq ~ up, then the rate A ~ 0, and there will be no further 
jumps in direction; the agent will continue traveling rightward forever. With these 
parameters, we will observe a front moving rightward. For example, for p = 0.017, 
q = 0.5, and u = 0.0165, this phenomenon will happen when the slope is larger 
than approximately 0.5, which is perfectly consistent with our simulations. On 
the other hand, for larger p, for example p = 1, this will not happen until the 
slope is very large (larger than about 30). In summary, for either faster dynamics 
{p larger) or smaller slopes (smaller p), we expect our diffusion approximation to 
be more accurate. This is consistent with the simulation results. 



-jj{z[t) +pq) , we have in the 
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